function [ct, lags,na, nb, dt] = autocorr(varargin)
%
% autocorr - calculate autocorrelations of input lists in excel file,
% assuming first column contains trace # and second contains event times.
% each trace is computed separately and all are combined to get an average
% cross correlation.
% P. Manis 10/25/2006 (in the air: Calgary to Chicago).
%



startdir = '/Users/pmanis/Desktop'; % define starting directory here
filefilter = '*.xls'; % for excel files


tmaxlag = 2000;

ct=[];
lags=[];
na = 0;
nb = 0;
dt = 0;

% read the excel file. First 2 columns contain information of interest.
thisdir = pwd;
cd(startdir);
[f p e] = uigetfile(filefilter, 'autocorr - find file');
if(isempty(f))
    cd(thisdir);
    return;
end;
fn = [p f e(1:end-1)];
if(~exist(fn))
    fprintf(1, 'File not found: \"%s\"\n', fn);
    cd(thisdir);
    return;
end;
type = xlsfinfo(fn);
if(~strmatch(type, 'Microsoft Excel Spreadsheet'))
    fprintf(1, 'File is not Excell spreadsheet\n');
    cd(thisdir);
    return;
end;
d = xlsread(fn);
cd(thisdir);
size(d)


tr  = d(:,1);
tru = unique(tr);
ntr = length(tru); % find out the unique event time lists.


maxlat = max(d(:,2)); % get max of events


minlat = 0;
c={}; lags = {};
dt = 1;
twin = tmaxlag;
maxlag = floor(tmaxlag/dt);

tmax = 0;
atmax = 0;
btmax = 0;
% set up 
for i = 1:ntr
    D1 = d(find(d(:,1) == tru(i)), 2);
    D2 = D1; % routine was designed for cross correlation - this makes it autocorr.
    t1 = max(D1/dt);
    t2 = max(D2/dt);
    if(t1 > atmax)
        atmax = t1;
    end;
    if(t2 > btmax)
        btmax = t2;
    end;
end;
if(atmax > btmax)
    tmax = atmax;
else
    tmax = btmax;
end;

zt = zeros(1, floor(tmax)+1); % build sample base - time .
ct = zeros(1, maxlag*2+1); % cumulative correlation result
na = 0; % number of points in first array
nb = 0; % and second (for cross correlation). 

% calculations start here...
for i = ntr
    D1 = d(find(d(:,1) == tru(i)), 2);
    D2 = D1; % routine was designed for cross correlation - this makes it autocorr.

    % first we convert the event lists to a list of 1's and 0's corresponding to the times
    % when the events occur...
    a = floor(D1/dt);
    a = a(find(a > minlat & a < maxlat));
    b = floor(D2/dt);
    b = b(find(b > minlat & b < maxlat));
    ia = zt;
    ia(a) = 1;
    ib = zt;
    ib(b) = 1;

    if(any(ia ~= 0) & any (ib ~= 0)) % make sure there is data to analyze
        na = na + length(find(ia == 1));
        nb = nb + length(find(ib == 1));
        [cb, lags] = xcorr(ia, ib, maxlag, 'none');
        ct = ct + cb; % get cumulative sum
    end;


end;
% ct = ct/size(D1,2); % average....
ct = 1000*ct/(dt*na); % normalizing by input event count ... to get rate out (in events per second).

% to look for periodicity, you need to delete the peak at 0 time...
ct(floor(end/2)+1) = 0;

% plot it.
h = findobj('tag', 'simcorrfig');
if(isempty(h))
    h = figure;
    set(h, 'tag', 'simcorrfig');
else
    figure(h);
    clf;
end;

bar(lags, ct);
set(gca, 'Xlim', [-tmaxlag tmaxlag]);

title(sprintf('Autocorrelation for %s (%s)',fn, date()), 'Interpreter', 'none');


